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Abstract: 


'  This  paper  defines  the  Difference  of  Low-Pass  (DOLP)  transform  and  describes  a  fast  algorithm  for 
its  computation.  The  DOLP  is  a  reversible  transform  which  converts  an  image  into  a  set  of  band-pass 
images.  A  DOLP  transform  is  shown  to  require  0(N*)  multiplies  and  produce  0(N  Log(N))  samples 
from  an  N  sample  image.  When  Gaussian  low-pass  filters  are  used,  the  result  is  a  set  of  images  which 

have  been  convolved  with  difference  of  Gaussian  ( DOG)  filters  from  an  exponential  set  of  sizes. 

>L-' 

A  fast  computation  technique  based  on  ^resampling"  is  described  and  shown  to  reduce  the  DOLP 
transform  complexity  to  0(N  Log(N))  multiplies  and  O(N)  storage  locations.  A  second  technique, 
"cascaded  convolution  with  expansion'",  is  then  defined  and  also  shown  to  reduce  the  computational 
cost  to  0(N  Log(N))  multiplies.  Combining  these  two  techniques  yields  an  algorithm  for  a  DOLP 
transform  that  requires  0(N)  storage  cells  and  requires  O(N)  multiplies. 

The  DOLP  transform  provides  a  basis  for  a  structural  description  of  gray-scale  shape.  Descriptions 
of  shape  in  this  representation  may  be  matched  efficiently  to  descriptions  of  shape  from  other  images 
to  determine  motion  or  stereo  correspondence.  Such  descriptions  may  also  be  matched  independent 
of  their  size  or  image  plane  orientation. 
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1  Introduction 


The  Difference  of  I.ow-Pass  (DOLP)  Transform  is  a  reversible  transform  which  converts  an  image  into  a 
set  of  band-pass  images.  These  band-pass  images  comprise  a  three  space  (the  DOI.P  space)  which  serves  as 
the  basis  for  an  efficient  technique  for  matching  descriptions  of  shape  [10]. 

The  band-pass  images  which  compose  the  DOI.P  space  arc  each  equivalent  to  a  convolution  of  the  image 
with  a  band-pass  filter,  bk.  Each  band-pass  filter  is  formed  by  a  difference  of  two  size-sealed  copies  of  a 
low-pass  filter,  gk/  and  gk. 

~  Sk-1  ~  Sk 

Each  low-pass  filter  gk  is  a  copy  of  the  low  pass  filter  gk  [  sealed  larger  in  size  by  a  scale  factor. 

In  the  following  sections  we  motivate  the  need  for  fast  computation  of  a  multi-resolution  description  of 
image  signals,  and  briefly  describe  a  representation  based  the  DOLP  transform.  This  representation  is  the 
topic  of  a  companion  paper  [11].  We  then  introduce  two  techniques  for  speeding  the  computation  of  a  DOLP 
transform.  A  fast  algorithm  based  on  these  techniques  is  described  below.  This  algorithm  reduces  the 
complexity  of  computing  a  DOLP  transform  from  0(N2)1  to  O(N)  multiplies  and  additions,  where  N  is  the 
number  of  sample  points  in  an  image. 


1.1  Motivation:The  Structural  Description  of  Images 

Interpreting  the  patterns  in  an  image  requires  some  form  of  matching.  If  the  interpretation  is  restricted  to 
two-dimensional  patterns,  this  matching  is  between  descriptions  of  shapes  in  the  image  and  object  models.  If 
the  interpretation  is  in  terms  of  three-dimensional  objects  then  techniques  for  matching  among  stereo  images 
or  motion  sequences  may  be  require  d,  in  cither  case,  the  matching  problem  is  simplified  if  descriptions  are 
compared  at  multiple  resolutions. 

Detecting  peaks  and  ridges  in  a  DOLP  Transform  provides  a  structural  description  of  the  grayscale  shapes 
in  an  image.  Matching  the  structural  descriptions  of  shapes  in  images  is  an  efficient  approach  to  determining 
the  three-dimensional  structure  of  objects  from  stereo  pairs  of  images  and  from  motion  sequences  of  images 
[13].  Matching  to  a  prototype  description  of  an  object  class  is  also  useful  for  recognizing  shapes  in  both 
two-dimensional  image  domains  and  three-dimensional  scene  domains  [3].  The  motivation  for  computing  a 
structural  description  is  to  spend  a  fixed  computational  cost  to  transform  the  information  in  each  image  into  a 
representation  in  which  searching  and  matching  arc  more  efficient.  In  many  eases  the  computation  involved 
in  constructing  a  structural  description  is  regular  and  local,  making  the  computation  amenable  to  fast 
implcmentauon  in  special  purpose  hardware. 

Several  researchers  have  shown  that  the  efficiency  of  searching  and  matching  processes  can  be  dramatically 
improved  by  performing  the  search  with  a  multi-resolution  hierarchy.  Moravee  [15]  has  demonstrated  a 
multi-resolution  correspondence  matching  algorithm  for  object  location  in  stereo  images.  Marr  and  Poggio 
[13]  have  demonstrated  correspondence  matching  using  edges  detected  by  a  difference  of  Gaussian  filters  at 


*The  symbol  0( )  is  pronounced  "order  of.  A  function.  g(n)  is  said  to  be  of  O(fln))  if  there  exists  a  constant,  c.  such  that  gfn)  <,  efin) 
for  ail  but  some  finite  (possible  empty)  set  of  nonnegative  values  for  n  [2]. 
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four  resolution.  Roscnfcld  and  Vandcrbmg[21]  have  described  a  two  stage  hierarchical  template  matching 
algorithm.  Hall  has  reported  using  a  multi-resolution  pyramid  to  dramatically  speed  up  correlation  of  aerial 
images  [12].  It  should  also  be  noted  that  Burt  has  recently  reported  using  cascaded  convolution  of  "Gaussian- 
Like"  kernels  to  construct  a  form  of  DOLP  transform  [4]. 

There  is  also  experimental  evidence  that  the  visual  systems  of  humans  and  other  mammals  separate  images 
into  a  set  of  "spatial  frequency”  channels  as  a  first  encoding  of  visual  information.  This  "multi-channel 
theory”  is  based  on  measurements  of  the  adaption  of  the  threshold  sensitivity  to  vertical  sinusoidal  functions 
of  various  spatial  frequencies  [7],  [22].  Adaption  to  a  sinusoid  of  a  particular  frequency  affects  only  the 
threshold  sensitivity  for  frequencies  within  one  octave.  ITiis  evidence  suggests  that  mammalian  visual  systems 
employ  a  set  of  band-pass  channels  with  a  band-width  of  about  one  octave.  Such  a  set  of  channels  would 
carry  information  from  different  resolutions  in  the  image.  These  studies,  and  physiological  experiments 
supporting  the  concept  of  parallel  spatial  frequency  analysis,  arc  reviewed  in  [6]  and  [23]. 


1 .2  The  Structural  Description  of  Shape  Based  on  the  DOLP  Transform 

The  DOLP  transform  provides  the  basis  for  a  representation  in  which  two-dimensional  gray  scale  shape  is 
described  by  a  tree  of  symbols  [10].  A  description  in  this  representation  contains  a  small  number  of  symbols 
at  the  root.  These  symbols  describe  the  global  (or  low-frequency)  structure  of  a  shape.  At  lower  levels,  this 
tree  contains  an  increasingly  larger  numbers  of  symbols  which  represent  more  local  events.  Finding  the 
correspondence  between  symbols  at  one  level  in  the  tree  constrains  the  possible  set  of  correspondences  at  the 
next  higher  resolution  level. 

The  description  is  created  by  detecting  local  positive  maxima  and  negative  minima  (peaks)  in  each  band¬ 
pass  image  of  a  DOLP  transform.  Local  peaks  in  the  DOLP  three-space  define  locations  and  sizes  at  which  a 
DOLP  band-pass  filter  best  fits  a  grayscale  pattern.  These  points  arc  encoded  as  symbols  which  serve  as 
landmarks  for  matching  the  information  in  images.  Peaks  of  the  same  sign  which  arc  in  adjacent  positions  in 
adjacent  band  pass  images  arc  linked  to  form  a  tree.  During  the  linking  process,  the  largest  peak  along  each 
branch  is  detected.  This  largest  peak  serves  as  a  landmark  which  marks  the  position  and  size  of  a  gray-scale 
form  (or  shape).  The  paths  of  the  other  peaks  which  arc  attached  to  such  landmarks  provide  a  further 
description  of  the  shape  of  the  form,  as  well  as  a  continuity  with  structural  forms  at  other  resolutions.  Further 
information  is  encoded  by  detecting  and  linking  two-dimensional  ridge  points  in  each  band-pass  image  and 
three-dimensional  ridge  points  within  the  DOLP  three-space. 


1 .3  A  Fast  DOLP  T ransform 

A  full  DOLP  transform  of  an  image  composed  of  N  samples,  produces  K  =  Logs(N)  band-pass  images 
composed  of  N  samples  each,  and  requires  0(N2)  multiplies  and  additions,  where,  S  is  a  "Scale  Factor"  which 
is  discussed  below.  Two  techniques  can  be  used  to  reduce  the  computational  complexity  of  die  DOLP 
transform:  "resampling"  and  "cascaded  convolution  with  expansion". 

Resampling  is  based  on  the  fact  that  the  filters  used  in  a  DOLP  transform  are  sealed  copies  of  a  band- 
limited  filter.  As  the  filter’s  impulse  response  becomes  larger  in  the  space  domain,  its  upper  cutofF  frequency 
decreases,  and  thus  its  output  can  be  resampled  with  coarser  spacing  widiout  loss  of  information.  The 
exponential  growth  in  the  number  of  filter  coefficients  which  results  from  the  exponential  scaling  of  size  is 


offset  by  an  exponential  growth  in  distance  between  points  at  which  the  convolution  is  computed.  The  result 
is  that  each  band-pass  image  may  be  computed  with  the  same  number  of  multiplications  and  additions. 
Resampling  each  band-pass  image  also  reduces  the  total  number  of  points  in  the  IX)l.P  space  from 
N  Logs(N)  samples  to  3N  samples. 

Cascaded  convolution  exploits  the  fact  that  the  convolution  of  a  Gaussian  function  with  itself  produces  a 
Gaussian  sealed  larger  in  sir.c  by  V2  .  This  method  also  employs  an  operation,  referred  to  as  "expansion",  in 
which  the  coefficients  of  a  filter  arc  mapped  into  a  larger  sample  grid,  thereby  expanding  the  size  of  the  filter. 
This  operation  can  be  used  without  introducing  distortion  under  certain  conditions  when  the  filter  is  band- 
limited.  and  is  to  be  convolved  with  a  band-limited  signal. 


1 .4  Organization  of  this  Paper 

Section  2  defines  the  DOLP  transform  and  shows  that  its  computation  requires  0(N2)  multiplies  and 
0(N  Log(N))  storage  locations.  Each  of  the  two  fast  computation  techniques  arc  described  and  their 
complexity  is  analyzed  in  section  3.  A  fast  algorithm  based  on  both  of  these  techniques  is  then  described  and 
shown  to  require  O(N)  multiplies  and  O(N)  Storage  locations.  An  example  is  then  presented  of  the  band-pass 
images  that  result  from  this  fast  algorithm  in  section  4. 

2  The  DOLP  T ransform 

This  section  defines  the  DOLP  transform  and  shows  that  its  computation  requires  0(N2)  multiplies  and 
0(N  l.og(N))  storage  locations.  This  is  followed  by  a  description  of  each  of  the  two  fast  computation 
techniques  and  an  analysis  of  the  computational  complexity  of  the  algorithms  based  on  each  technique.  A 
fast  algorithm  based  on  both  of  these  techniques  is  then  described  and  shown  to  require  O(N)  multiplies  and 
0(N)  Storage  locations. 


2.1  The  DOLP  Transform  Definition 

The  DOLP  transform  expands  an  N  x  N  image  signal  p(x,y)  into  Logs(N)  band-pass  images  %^x,y).  Each 
band-pass  image  is  equivalent  to  a  convolution2  of  the  image  p(x,y)  with  a  band-pass  impulse  response, 
b/x.y). 

%Jx,y)  =  p(x,y)  *  bjfx.y)  (1) 

The  DOLP  transform  is  illustrated  in  the  data  flow  graph  shown  in  figure  1. 

For  k =0.  the  band-pass  filter  is  formed  by  subtracting  a  circularly  symmetric  low-pass  filter  gjx.y)  from  a 
unit  sample  positioned  over  the  center  coefficient  at  the  point  (0,0). 

bjx,y)  =  S(x,y)  -  g0(x,y)  (2) 


c  filters  described  in  this  paper  are  all  non- recursive  finite  impulse  response  (FIR)  filters.  Convolutions  are  computed  for  each 
ample  point  in  the  image:  when  the  filter  coefficients  extends  beyond  the  tdge  of  the  image,  a  default  border  value  (typically  0  )  is 
supplied  in  place  of  the  image  value 


b0  (x,  l ) 
bj(x,  y ) 
ba(x,  y ) 


Figure  1:  The  Difference  of  Low-Pass  (DOLP)  Transform 

This  data  flow  graph  illustrates  the  computational  process  for  a  DOLP  transform.  The 
transform  produces  Logs(N)  band-pass  images.  Each  band-pass  image  is  produced  by 
convolving  the  image  with  a  band-pass  impulse  response  (filter)  which  is  a  size-scaled 
copy  of  a  prototype  filter.  This  prototype  is  formed  from  a  difference  of  two  size-scaled 
copies  of  a  low-pass  filter. 

The  filter  b0(x,y)  gives  a  high-pass  image,  ?S>0(x,y).  This  image  is  equivalent  to  the  result  produced  by  the 
edge  detection  technique  known  as  "unsharp  masking"  [20], 

(x,y)  =  P(x,y)  *  (  fi(x,y)  -  ga(x,y)  )  (3) 

=  P(x,y)  -  (p(x,y)  *  g0(x,y)) 


For  band-pass  levels  1  ^  k  <  K  the  band-pass  filter  is  formed  as  a  difference  of  two  size-scaled  copies  of  the 
low-pass  filter. 

bjx,y)  =  gh/x,y)  -  g/x,y)  (4) 

Each  low-pass  filter,  gk(x,y)  is  a  copy  of  the  circularly  symmetric  low-pass  filter  g0(x,y)  scaled  larger  in  size 
by  a  factor  raised  to  the  k111  power.  Thus  for  each  k,  the  band-pass  impulse  response,  bk(x.y),  is  a  size  scaled 
copy  of  the  band-pass  impulse  response,  bk/(x,y).  This  property  is  necessary  for  the  configuration  of  peaks  in 
a  DOLP  transform  of  a  shape  to  be  invariant  to  the  size  of  the  shape  [10]. 


The  scale  factor  is  an  important  parameter  which  affects  several  aspects  of  the  DOLP  transform.  For  a 


two-dimensional  DOLP  transform,  this  scale  factor,  denoted  S2,  has  a  typical  value  of  V2 .  In  the  ease  of  a 
one-dimensional  1)01  ,P  transform,  the  scale  factor  is  denoted  and  has  a  typical  value  of  2.  This  scale 
factor  is  discussed  again  at  the  end  of  this  section. 


For  two-dimensional  circularly  symmetric  filters  which  are  defined  by  sampling  a  continuous  function,  size 
scaling  can  defined  as  increasing  the  density  of  sample  points  over  a  fixed  domain  of  the  function.  In  the 
«  Gaussian  filter,  this  has  the  effect  of  increasing  the  standard  deviation,  a,  relative  to  the  image  sample  rate. 

In  principle  the  DOLP  transform  can  be  defined  for  any  number  of  band-pass  levels  K.  A  convenient  value 
ofKis 

K  =  I.ogs(N) 

where  S  is  equal  to  the  sample  distance  Sx  for  a  one-dimensional  DOLP  transform,  or  the  square  of  the 
sample  distance  S2  for  a  two-dimensional  DOLP  transform. 


This  value  of  K  is  the  number  of  band-pass  images  that  result  if  each  band-pass  image,  ®k,  is  resampled  at 
a  sampling  distance  of  S2.  With  this  resampling,  the  K*  image  contains  only  one  sample. 

The  DOLP  transform  is  reversible.  The  original  image  may  be  recovered  by  adding  all  of  the  band-pass 
images,  plus  a  low-pass  residue.  This  low  pass  residue,  which  has  not  been  found  to  be  useful  for  describing 
the  image,  is  obtained  by  convolving  the  lowest  frequency  (largest)  low-pass  filter,  g^fx.y)  with  the  image. 

K-l 

p(x,y)  =  {p(x,y)  *  g/x,y))  +  2Z  ^k(x,y)  .  (5) 

k=0 


Reversibility  proves  that  no  information  is  lost  by  the  DOLP  Transform. 

Because  convolution  and  subtraction  are  associative  the  DOLP  transform  may  also  be  computed  by 
convolving  the  original  image  with  an  exponentially  size-sealed  set  of  low-pass  filters  and  subtracting  each 
low-pass  image  from  the  next  to  form  the  set  of  band-pass  images.  This  technique  is  illustrated  in  figure  2. 
One  of  the  fast  computational  techniques  for  a  DOLP  transform  which  are  described  below  is  based  on  the 
idea  of  computing  the  convolutions  of  the  image  with  progressively  larger  low-pass  filters  which  are 
implemented  as  a  cascade  of  convolutions  with  small  low-pass  filters. 

2.2  Discussion:  The  Scale  Factor 

The  parameter  S,  used  to  determine  the  number  of  levels,  is  crucial  to  both  the  scaling  of  low-pass  filters 
and  resampling  of  the  band-pass  and  low-pass  images.  These  two  ideas  arc  related  when  peaks  and  ridges 
from  the  DOLP  transform  are  to  be  used  to  describe  the  shape  of  a  form  so  that  it  can  matched  independent 
of  the  size  of  the  form.  In  such  an  application  it  is  important  that  the  density  of  samples  be  a  fixed  fraction  of 
the  size  of  the  band-pass  impulse  response.  It  is  convenient  to  define  a  single  variable,  S  =  S2  =  to 
simplify  the  expression  for  K  as  well  as  for  some  of  the  analysis  given  below. 
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Figure  2:  The  Difference  Method  for  Computing 
the  Difference  of  Low-Pass  (DOLP)  Transform 

Because  convolution  and  subtraction  are  associative  the  DOLP  transform  may  also  be 
computed  by  convolving  the  original  image  with  an  exponentially  size-scaled  set  of 
low-pass  filters  and  subtracting  each  low-pass  image  from  the  next  to  form  the  set  of 
band-pass  images.  The  data  flow  graph  for  this  process  shows  the  reversibility  of  the 
DOLP  transform.  This  approach  is  also  the  basis  for  a  fast  computation  technique  for 
the  DOLP  transform  called  "Cascaded  filtering  with  expansion".  With  this  technique 
the  sequence  of  low-pass  images  are  obtained  by  repeated  convolution  with  a  small 
kernel  filter. 

Marr  [14]  argues  that  a  value  of  S2=  1.6  js  "optimum"3  for  a  difference  of  Gaussian  band-pass  filter.  For 
two-dimensional  signals  the  value  S2=  V2  has  virtually  the  same  effect  while  providing  some  additional 
benefits. 


Marr  calls  this  value  optimum  in  the  sense  that  it  simultaneously  minimizes  S2  while  maximizing  the  energy  in  the  filter.  A  curve  of 
filter  energy  with  respect  to  ratio  of  standard  deviations  exhibits  a  "knee”  in  the  region  of  1.6.  [14].  For  smaller  ratios  the  energy  of  the 
resulting  filter  falls  rapidly,  while  for  larger  values  it  is  nearly  constant 
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The  most  important  benefit  of  using  S2=  V2  is  that  V2  is  the  smallest  naturally  occuring  resample 
distance  on  a  two-dimensional  cartesian  grid.  Thus  by  using  S2  =  .\/2  we  can  resample  each  band-pass  image 
at  a  distance  that  is  a  constant  fraction  of  the  band-pass  filter  size.  This  yields  a  configuration  of  peaks  and 
ridges  in  a  DOl.P  transform  that  is  invariant  to  the  size  of  a  shape,  except  for  cyclic  distortions  due  to 
sampling  effects.  Such  descriptions  of  shapes  can  be  matched  independent  of  the  size  of  the  shape. 

An  additional  benefit  from  using  S2=  V2  comes  from  the  Gaussian  auto-convolution  scaling  property. 
When  a  Gaussian  function  is  convolved  with  itself  the  result  is  the  Gaussian  function  sealed  larger  in  size  by 
V2  .  We  will  show  below  that  this  property  can  be  used  to  greatly  reduce  the  computational  cost  of  a  DOLP 
transform. 


2.3  Complexity  of  DOLP  Transform 

In  this  section  we  derive  formulae  for  the  memory  requirements  and  computational  costs  of  the  DOLP 
transform.  A  first  step  in  obtaining  these  quantities  is  the  calculation  of  the  number  of  coefficients  in  each 
filter.  We  do  this  for  both  the  one  and  two-dimensional  eases  and  then  produce  a  general  result  that  holds  in 
both  cases. 

Let  Rk  refer  to  the  radius  of  the  filter,  and  let  Xk  refer  to  the  number  of  coefficients,  for  both  the  one  and 
two-dimensional  cases.  Also,  let  Sj  refer  to  the  size  scaling  factor  for  the  one-dimensional  filters  and  S2  refer 
to  this  factor  for  two-dimensional  filters,  as  above. 

In  the  one-dimensional  case,  the  number  of  coefficients  is  specified  by  the  radius  of  the  filter. 

Xk  =  2  Rk  +  1 

The  radii  at  each  band-pass  level  k  are  related  to  the  radius  R0  of  the  smallest  level  by 

Rk  =  Ro 

Thus  the  number  of  coefficients  for  the  k*  band-pass  filter  is 

Xk  =  (X0-l)S{  +  1 

Since  X0  »  1  we  can  simplify  the  mathematics  by  using  the  approximation: 

Xk  ~  X0  s2 

In  the  case  of  the  two-dimensional  filters  for  images,  the  low-pass  filter,  g0(x,  y),  is  specified  to  be  circularly 
symmetric.  If  the  coefficients  arc  nonzero  for  all  points  (x,  y)  such  that  x2  +  y2  <  R2  then, 

X0  a  7T  R2 

This  approximation  becomes  more  accurate  as  R0  increases. 
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The  band-pass  filters  for  levels  1  through  K-l  are  specified  to  be  size  sealed  copies  of  the  level  0  filter. 
Each  filter  is  to  be  sealed  larger  in  size  by  a  factor  of  S2.  Thus  Rk  is  related  to  R0  by 
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Rk  =  R„S![ 

As  a  result  the  number  of  coefficients  at  level  k  is 
Xk  s Xa Sf 

If  we  define  the  variable,  S,  such  that  S  =  Sj  =  as  before,  then  in  both  the  one-dimensional  and 
two-dimensional  ease, 

Xk  s  X0  Sk 

This  approximation  becomes  more  accurate  as  k  increases. 

As  described  above,  the  DOLP  transform  is  defined  to  produce  band-pass  levels  0  through  K-l,  where  K 
is 

K  =  Logs(N) 

Since  the  DOLP  transform  produces  K  band-pass  images  of  N  samples  each,  the  memory  requirement  M  is 

(6) 

M  =  N  K  =  N  LoggCN) 

The  number  of  multiplies  for  producing  each  band-pass  image  is  proportional  to  the  number  of  coefficients 
in  the  filter  for  that  level.  rlhe  total  number  of  multiplies  for  the  convolutions,  denoted  C  (for  cost),  is  given 
by: 

C  s  N  X„(l  +  1  +  S  +  S2  +.  +  SK_1) 


*X0N(1  +  D  Sk) 


«  X9  N  ( 1  + 


Using  our  typical  value  S  =  2, 


i  -,K 

1  + - =  2k 


and  the  cost  becomes: 


C  s  XB  N  2k  =  X.  N 


and  thus 


C  ss  X-  N2 


Figure  3: 

Techniques  for  Reducing  the  Cost  of  a  DOLP  Transform 

Two  independent  techniques  can  be  used  to  reduce  the  computational  cost  of  a 
DOLP  transform:  Resampling  and  Cascaded  Convolution  with  Expansion.  These  two 
techniques  can  be  combined  to  produce  an  algorithm  which  for  computing  a  DOLP 
transform  in  0(N)  multiplies  which  requires  O(N)  storage  cells. 


3  Fast  Computation  Techniques 


Wc  have  developed  two  independent  techniques  to  reduce  the  computational  cost  of  a  1)01. P  Transform. 
Kadi  of  these  techniques  reduces  the  cumber  of  multiplies  and  additions  for  an  N  sample  1)01. P  transform 
from  0(N2)  multiplies  to  0(N  I.og(N)  )  multiplies  and  additions.  Combined,  these  techniques  allow  the 
DOI.P  transform  to  be  computed  with  0(N)  multiplies  and  0(N)  additions. 

The  two  techniques  are: 

•  Resampling:  Computation  of  the  band-pass  images  at  resample  points  which  arc  spaced  at  a  fixed 
fraction  of  the  filter  radius. 

•  Cascade  Convolution  with  Expansion:  Use  of  the  autoconvolution  scaling  property  of  the 
Gaussian  low-pass  filter  and  a  remapping  of  the  filter  coefficients  to  obtain  the  impulse  response 
of  a  larger  filter  from  a  cascade  of  small  filters. 

These  two  techniques  may  be  applied  independently  to  reduce  the  computational  cost  of  a  DOI.P  transform, 
as  illustrated  in  figure  3.  When  combined,  these  two  techniques  provide  an  algorithm  which  will  compute  a 
DOI.P  transform  in  0(N)  multiplies  with  a  storage  requirement  of  0(N)  cells.  In  the  following  sections  we 
describe  algorithms  for  computing  a  DOLP  transform  based  on  each  of  these  techniques  separately.  We  then 
describe  the  algorithm  which  employs  both  techniques. 

This  section  begins  with  a  discussion  of  resampling  a  cartesian  two-dimensional  signal  at  a  distance  of  V2 . 
A  linear  systems  modcl_  for  such  resampling  is  presented.  We  then  describe  the  Sampled  DOLP  transform, 
and  show  that  with  \Jl  resampling,  a  DOLP  transform  can  be  computed  with  0(N  Log(N))  multiplies  and 
that  this  DOLP  transform  can  be  stored  in  0(N)  storage  cells. 

.  We  then  discuss  the  scaling  property  of  the  Gaussian  filter,  and  show  that  a  Gaussian  impulse  response  of 
size  SV2  can  be  formed  by  convolving  a  Gaussian  filter  of  size  S  with  itself.  This  technique  is  referred  to  as 
cascaded  convolution.  A  second  scaling  operation  known  as  the  expansion  operator  is  then  introduced.  We 
show  that  a  combination  of  expansion  and  cascaded  convolution  can  also  be  used  to  compute  a  DOLP 
transform  of  an  N  sample  image  in  0(N  Log(N))  multiplies. 

Finally,  these  two  techniques  are  combined  to  produce  an  algorithm  which  will  compute  a  DOLP 
transform  which  requires  0(N)  samples  in  0(N)  multiplies.  This  technique  is  referred  to  as  "Cascaded 
Convolution  with  Expansion  and  Resampling. 


3.1  Resampling 

The  number  of  samples  that  is  needed  to  represent  a  discrete  signal  is  determined  by  the  frequency  content 
of  that  signal.  As  Nyquist  demonstrated,  [16],  a  signal  which  has  been  convolved  with  a  filter  which 
attenuates  the  higher  frequency  components  may  be  represented  by  a  smaller  number  of  samples.  Very  little 
information  is  lost  when  a  band-limited  signal  isjresamplcd  because  the  original  samples  may  be  recovered  by 
interpolation.  In  this  section  we  describe  the  V2  sampling  operation  and  then  present  the  algorithm  for  the 
sampled  DOLP  transform. 
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3.1.1  Sampling  at  \/2 
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Figure  4:  Example  of  S^j{p(x,y)}  (  Circles) 
andS^/jtS^/jtpfx.y)}}  (Squares) 

Applied  to  a  Cartesian  Sample  Grid  (  Dots  ) 

The  smallest  distance  between  sample  points  on  a  two-dimensional  cartesian  grid  which  is  larger  than  1  is 
V2  ,  the  distance  between  diagonally  adjacent  elements.  A  two-dimensional  signal  may  be  resampled  at  this 
sample  distance  by  removing  every  other  diagonal,  as  illustrated  by  the  circles  in  Figure  4.  We  refer  to  this 
process  as  V2  sampling,  denoted  S^/jt}.  V2  sampling  reduces  the  number  of  sample  points  in  a  two- 
dimensional  signal  by  1/2.  A  second  application  of  V2  'sampling  will  produce  a  two-dimensional  signal 
which  has  samples  spaced  at  a  distance  of  2  on  the  original  grid,  as  shown  by  the  boxes  in  figure  4. 

The  points  on  the  V2  sample  grid  may_be  detected  by  a  simple  test  using  the  modulus  function  ( denoted 
here  as  "mod" ).  Sample  points  on  the  V2  grid  are  those  points,  (x,y)  which  satisfy  the  relationship 

x  mod  2  =  y  mod  2. 

This  is  the  sample  function  applied  to  level  2  of  the  sampled  DOLP  transform.  A  second  application  of 
V2  sampling  produces  a  sample  grid  with  a  minimum  distance  of  2  between  samples.  These  points  are  those 
for  which 

x  mod  2  =  0  and  y  mod  2  =  0 
This  is  the  sample  grid  for  level  3  of  the  sampled  DOLP  transform. 

In  general,  each  level  k,  for  2  <_k  <  K-l,  of  the  sampled  DOLP  transform  will  have  a  sample  grid 
produced  by  k-l  applications  of  V2  sampling.  Those  levels  for  which  k  is  even  will  have  sample  points 
defined  by 

x  mod  2^'^/2  =  y  mod  2*k'^/2 

Those  levels  for  which  k  is  odd  will  have  points  which  are  given  by 
x  mod  2(k'l)/2  -  0 
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and 

y  mod2(k’1^2  =  0 

Such  sampling  may  be  done  for  any  value  S  which  is  a  distance  between  points  on  the  original  sample  grid. 
For  example,  if  we  select  points  that  are  separated  by  a  distance  in  the  x  dimension  of  2  and  in  the  y 
dimension  of_l.  then  our  resample  distance  is  S2  =  V2Z+1  =  V5  .  If  a  two-dimensional  scale  factor  other 
than  S2  =  \/2  is  used,  the  value  S  =  S2  must  be  substituted  for  the  2  in  the  above  expressions.  In  this  ease 
the  size  of  the  low-pass  filters  should  be  sealed  by  this  same  factor  if  the  DOI.P  transform  is  to  be  used  to 
produce  a  description  of  shape  that  can  be  matched  at  any  size. 

3.1.2  Linear  Systems  Model  for  Resampling 


Figure  5:  Nyquist  Boundaries  for  Successive  applications  of  V2  Resampling 

The  effects  of  resampling  are  best  described  in  the  spatial  frequency  domain.  Let  us  describe  the  transfer 
function  (  Discrete  Time  Fourier  Transform )  of  a  two  dimensional  function,  h(x,y)  as  [17] 

00  oo 

H(u,v)  =  £  £  h(x,y)e-juxe*jvy  (8) 

u  =-00  V  =  -00 

The  continuous  variables  u  and  v  are  referred  to  as  the  spatial  frequency  variables.  Figure  5  shows  the  the 
range  of  unique  spatial  frequency  components  in  the  (u,  v)  plane  that  is  generated  by  the  transfer  function  of 
a  two-dimensional  signal.  A  two-dimensional  function  sampled  on  a  cartesian  grid  has  a  transfer  function 
which  is  unique  within  the  square  region  of  the  (u,  v)  plane  bounded  by  (  ±ir,  ±ir  ).  The  boundaries  of  this 
region  are  referred  to  as  the  Nyquist  boundaries.  The  resampling  operation  Sy^-{.}  generates  a  new  Nyquist 
boundary,  shown  by  the  diamond  shaped  region  in  Figure  5.  The  V2  resampling  operation  has  the  effect  of 
"folding’’  or  aliasing  any  signal  energy  outside  this  new  Nyquist  boundary.  This  folded  signal  energy  is  added 
to  the  signal,  and  appears  as  energy  at  a  lower  frequency.  Such  a  distortion  is  not  reversible  and  will 
introduce  errors  when  used  with  techniques  which  are  based  on  detecting  peaks  and  ridges. 


Aliasing  is  minimized  by  filtering  the  two-dimensional  sequence  so  that  there  is  very  little  signal  energy 
outside  the  Nyquisc  boundary  when  the  signal  is  resampled.4  'I"his  minimizes  the  reflected  signal  energy  that 
results  in  aliasing.  Mathematically,  the  operation  is  modelled  as  first  convolving  the  signal  with  a  band- 
limited  filter,  and  then  selecting  only  the  subset  of  points  at  which  the  filter  signal  is  resampled.  For 
implementation  on  a  serial  processor,  the  computational  cost  may  be  reduced  by  only  evaluating  the 
convolution  expression  at  those  points  where  the  filter  is  centered  over  the  resample  points.  This  "resampled 
convolution"  is  illustrated  by  the  function  Syy{}  placed  in  boxes  adjacent  to  the  convolution  boxes  in  Figure 
*  6. 

3.1.3  Complexity  of  the  Sampled  DOLP  Transform 

A  convolution  may  be  expressed  as  a  sequence  of  inner  products  of  the  filter  coefficients  with 
neighborhoods  of  the  signal.  By  only  computing  these  inner  products  for  the  instances  where  the  filter  is 
centered  over  resampled  points,  it  is  possible  to  reduce  the  computational  complexity  of  a  DOl.P  transform  to 
0(N  I.og(N) ).  In  such  a  Sampled  DOl.P  transform,  the  distance  between  resample  points  increase  by  the 
same  scale  factor  as  the  band-pass  filters.  The  computational  complexity  and  memory  requirements  for  the 
Sampled  DOLP  Transform  may  be  evaluated  by  considering  the  steps  in  the  algorithm.  In  this  section  we 
present  such  an  analysis  for  any  value  of  S. 

The  band-pass  signals.  S0(x,y)  and  Sjfx.y),  arc  computed  as  described  for  the  DOLP  transform,  requiring 
X„N  and  S  X0  N  multiplies  respectively.  <3&2(x,y>  is  computed  only  for  sample  points  in  p(x,y)  on  alternate 
diagonals.  The  convolution  at  level  2  is  with  a  filter  with  X0  S2  coefficients.  However,  the  convolution  is  only 
evaluated  at  the  N/S  sample  points  on  alternate  diagonals.  Thus  the  cost  is  S  X0  N  multiplies,  as  it  was  with 
level  1.  At  level  3.  the  band-pass  impulse  response  is  computed  for  sample  points  spaced  at  a  distance  of  Sj. 
There  arc  N/S2  such  points  and  the  filter  has  X0  S3,  so  the  cost  is  S  X„  N  multiplies. 

In  general  at  each  level  k  ,  for  2  <  k  <  K-I,  the  band-pass  filter  has  X0  Sk  coefficients,  and  the 

convolution  is  computed  at  N/S*kl*  sample  points,  for  a  cost  of  X0  S  N  multiplies  and  additions  at  each 
band-pass  level.  Since  there  are  K  =  Logs(N)  band-pass  levels,  the  total  cost  is 

C  =  X0  N  ( S  (  Logs(N)  -  1 )  + 1 )  multiplies  and  additions  (9) 

Band-pass  levels  0  and  1  each  have  N  samples.  For  levels  2  through  K-l  the  number  of  memory  cells 

required  drops  by  a  factor  of  S  for  each  level. 

M  =  N  (1  +  1  +  1/S  +  1/S2  +  1/S3  +  ...) 

K-l 

=  N(i  +  Ei 

k=0s 


4It  is  impossible  to  filter  a  sequence  with  a  finite  duration  filter  so  that  a  frequency  region  of  any  finite  size  is  identically  zero  [18]. 
However,  a*signal  can  be  filtered  so  that  there  is  an  arbitrarily  small  response  to  a  range  of  frequencies. 


b0(x>  tfl 

bj_(x,  y) 

b0(x,-y) 

S^(.} 

b3(x,  jf ) 

bH(x,  y ) 

S^3(.) 

Figure  6:  Data  Flow  Graph  for  Algorithm  for  Computing  Resampled  DOLP  Transform 

The  boxes  marked  with  Sy^fkjf.]  following  each  convolution  indicate  that  the 
convolution  is  computed  only  for  resample  points  specified  by  V2  resampling  at  level 
k.  ( See  text ) 


=  N(1+ - -) 


For  our  typical  value  of  S  =  2, 

M  =  N  (1  +  1  +  1/2  +  1/4  +  1/8  +...) 
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3.2  Cascaded  Convolution  with  Expansion 

Much  of  the  cost  of  a  IXDI-P  transform  results  from  the  large  number  of  coefficients  in  tlic  filters  for  larger 
values  of  k.  Resampling  compensates  for  the  exponential  increase  in  the  filter  size  by  an  exponential  increase 
in  the  space  between  sample  points.  A  second  technique  for  reducing  the  complexity  of  a  DOI.P  transform  to 
0{N  l.og  N)  multiplies  is  referred  to  as  "cascaded  convolution  with  expansion".  Phis  method  exploits  two 
mathematical  properties:  (1)  the  size-sealed  replication  of  the  Gaussian  functional  form  as  the  result  of  the 
convolution  of  a  Gaussian  function  with  itself.  (2)  a  scaling  operation  that  is  based  on  remapping  the 
coefficients  of  a  filter  into  a  new  sample  grid,  leaving  zero  or  undefined  samples  between  the  samples  of  the 
remapped  filer. 

In  the  following  sections  we  first  discuss  the  two-dimensional  circularly  symmetric  Gaussian  filter,  and  its 
properties  under  convolution.  We  then  describe  the  expansion  operator  and  the  algorithm  for  cascaded 
convolution  with  expansion,  together  with  an  analysis  of  its  complexity. 

3.2.1  The  Two-Dimensional  Circularly  Symmetric  Gaussian  Filter 

In  cascaded  convolution,  an  impulse  response  covering  a  large  support  is  obtained  by  repeatedly 
convolving  the  signal  with  copies  of  an  impulse  response  over  a  smaller  support.  This  algorithm  will  only 
produce  size  sealed  copies  of  the  low-pass  impulse  response  if  Gaussian  low-pass  filters  arc  used.  This  may  be 
shown  by  the  Gaussian  autoconvolution  scaling  property,  described  below. 

The  Gaussian  function  is  most  commonly  known  in  its  one-dimensional  form 

g(t)  A  1  Q-(t-ii)2/la2 

aV2ir 

where  fi  is  refered  to  as  the  mean  and  a  as  the  standard  deviation. 

The  term  l/a s/2tt  scales  the  Gaussian  function  so  that  it  has  unit  area. 

A  discrete  two-dimensional  Gaussian  filter  may  be  obtained  by  assuming  a  zero  mean  and  applying  the 
substitution 


r2 

a2  = — ,  and 
2  a 

The  coefficients  are  then  obtained  by  sampling  the  continuous  function  at  the  points  given  by  the  discrete 
variables  x  and  y  where  t2  =  x2  +  y2  <  R2. 

Implicit  in  this  filter  is  multiplication  by  a  uniform  circular  window  (or  aperture  or  support),  the  disk 

cR(x,y)  A  |  1  for  x2+y2<R2 
\  0  otherwise 
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To  control  the  filter  gain,  the  filter  coefficients  arc  normalized  so  that  they  sum  to  1.0.  This  is  done  by 
summing  the  coefficients  and  then  dividing  each  coefficient  by  the  sum. 

Thus  the  normalized  two-dimensional  Gaussian  low-pass  filter  defined  over  a  circular  support  is  given  by: 
g0(x.y)  =  (i/A)cR(x,y>e'a(x  +y2)/R 

Where  A  is  a  gain  factor  given  by 

A  =  X-  X]  cR(x,y)  e'a(x  +y2)/R 
|x|<R  |y|<R 

The  circularly  symmetric  function  cR(x,y)  has  a  transfer  function  [19] 


CR(u,v)  = 


2ffRJ,(RVu2+vz) 

Vu2+v2 


where  Jj()  is  the  first  order  Bessel  function. 


The  Gaussian  filter  ga(x,y)  has  a  transfer  function  which  is  a  Gaussian  function  convolved  with  the  transfer 
function  of  its  aperture  (or  support)  [19]. 


G.(u.»)  =  — cR(u,v) .  <-^)  e  -*V+ A«. 

A  RVir 

An  experimental  procedure  has  shown  that  the  parameters  R  =  4.0  and  a  =  4.0  work  quite  well  for 
cascaded  convolution  with  expansion  [10].  With  these  parameters,  the  transfer  function  of  the  impulse 
response  has  its  first  zero  crossing  in  a  circle  of  radius  approximately  equal  to  ir.  This  gives  a  filter  with  a 
pass-band  and  transition  region  which  just  fits  within  the  Nyquist  boundary. 

The  Gaussian  is  the  only  two-dimensional  function  which  is  both  circularly  symmetric  and  separable  into 
one-dimensional  components.  If  the  Gaussian  kernel  is  multiplied  by  a  square  support  rather  than  a  circular 
disc,  then  the  entire  impulse  response  can  be  separated  into  a  cascade  of  one-dimensional  components.  In 
this  case,  the  correlation  operation  can  be  implemented  with  significantly  fewer  multiplies  by  replacing  the 
convolution  with  a  (2R  +  1)  x  (2R  +  1)  circular  filter  by  two  convolutions  with  2R+ 1  point  one-dimensional 
filters  (  one  for  each  dimension).  This  requires  a  total  of  only  4R  +  2  multiplications  for  each  picture  point 
instead  of  4R2+4R  +  l  multiplications  [17].  The  square  support  degrades  the  circular  symmetry  of  the 
function.  The  result  is  some  additional  aliassing  along  the  axes  when  the  filtered  sequence  is  resampled. 


3.2.2  Cascaded  Convolution 

It  can  be  easily  shown  that  a  Gaussian  function  convolved  with  itself  yields  a  Gaussian  function  whose 
standard  deviation  is  V2  larger  than  the  original  function.  For  example,  in  one  dimension,  the  convolution 
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1  g-t2/2q2  * _ - _ g 

aV  lit  <jy/2w 


1  -tW 


may  also  be  expressed  as  the  product  of  Fourier  transforms 

2  2.-  2  2..  2  2 
Q-a  a  /2  9  g-cr  <J  /2  _  g-«  <0 


#  The  inverse  Fourier  transform  of  this  product  is 


alVir 

Returning  to  standard  form  requires  the  substitution 
<Tj  =  2a2  or  ff1  =  V2  <r. 

Thus  the  standard  deviation,  and  hence  the  function  width,  have  been  expanded  by  a  factor  of  V2 .  Note 
also  thauiutocon volution  preserves  the  unit  area  normalization;  the  amplitude  has  been  multiplied  by  a  factor 
of  1/V2  .  The  discrete  Gaussian  filter  is  of  finite  extent,  and  thus  is  not  an  exact  Gaussian.  For  this  reason 
the  Gaussian  scaling  property  only  holds  as  an  approximation  for  the  discrete  Gaussian  filter. 

Cascaded  convolution  provides  an  inexpensive  method  to  obtain  the  convolution  of  an  image  with  g^x.y). 
That  is,  low-pass  image  1  is  obtained  from  low-pass  image  0  by  a  second  convolution  with  g0(x,y),  yielding  the 
effective  filter, 

g,(x.y)  =  g0(x,y)  *  g0(x,y) 

However,  low-pass  image  2  then  requires  two  additional  convolutions  with  g0(x,y ).  and  low-pass  image  3 
requires  four  more  such  convolutions  with  g0(x,y).  This  exponential  growth  may  be  averted  by  resampling 
each  low-pass  image  by  V2  before  the  next  convolution,  or  by  expanding  the  ga(x,y)  onto  a  larger  sample 
grid  with  the  V2  expansion  operator. 

3.2.3  The  Expansion  Operator 

In  addition  to  cascaded  convolution  we  also  employ  a  technique  refered  to  as  "expansion"  in  the  algorithm 
described  below.  Expansion  is  possible  because  we  are  using  low-pass  filters  that  are  designed  with  a  high- 
frequency  stop  band.  These  filters  attenuate  the  spurious  high-frequency  signals  created  by  the  "expansion" 
operation. 

The  expansion  operation  is  a  spatial  remapping  of  the  samples  of  a  filter  so  that  the  distance  between 
samples  is  altered.  This  remapping  docs  not  affect  the  number  of  samples  in  a  filter  or  the  values  of  these 
samples.  Algorithm  arc  described  below  in  which  expansion  is  used  as  a  method  of  scaling  the  impulse 
response  larger  in  size  by  a  factor  of  Expansion  by  V2  is  necessary  in  order  to  convolve  a  filter  with  an 
image  which  has  been  resampled  to  a  V"2  sample  grid,  as  is  required  when  cascaded  convolution  is  used  with 
\Jl  resampling.  However,  it  is  also  possible  to  use  this  expansion  to  size-scale  a  filter  which  is  to  be 
convolved  with  a  conventional  cartesian  grid.  The  only  restriction  is  that  the  high  frequency  energy  generated 
by  expansion  must  be  attenuated  by  other  filters  in  the  cascade. 
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Hie  expansion  operation  may  be  modeled  as  a  spatial  scaling  followed  by  a  resampling.  A  simpler  analysis 
can  be  performed  by  considering  the  spacing  between  coefficients.  Both  analysis  produce  the  same  result: 
flic  transfer  function  of  the  filter  is  sealed  smaller  in  frequency  by  the  expansion,  and  copies  of  the  transfer 
function  appear  reflected  over  a  new  Nyquist  boundary  imposed  by  the  space  between  samples.  The 
conditions  under  which  expansion  can  be  used  without  distorting  the  image  arc  always  the  same.  The 
composite  cascade  filter  must  have  a  very  high  attenuation  everywhere  outside  of  the  new  Nyquist  boundary 
of  die  sample  grid  onto  which  the  filter  coefficients  arc  mapped. 

Let  us  define  (x.y)  as  points  in  the  cartesian  grid  in  which  a  filter  is  defined,  and  (xc,yc)  as  the 
corresponding  points  in  a  V2  grid  onto  which  the  filter  is  remapped.  A  single  application  of  the 
V2  expansion  operation  maps  each  row  from  a  filter  on  a  cartesian  sample  grid  into  every  other  diagonal  of 
the  V2  grid.  This  mapping  takes  each  coefficient  from  point  (x.y)  of  a  filter  g(x,y)  and  places  it  at  point 
(x-y,x  +  y)  of  a  filter  ge(xe,yc).  Points  of  gc(xc,ye)  which  receive  no  coefficient  under  this  mapping  arc  declared 
to  be  undefined  or  zero. 

Let  us  define  this  mapping  as  the  function  }•  Since 

xe  =  x  -  y 
ye  =  x  +  y 

we  obtain 


xe+ye 

x  =  -£ — S 
2 

and 


-xe+ye 

y  = _ _ s 


Thus  this  function  may  be  defined  by 

E x/2  {g(x.y)}  -  ge(xe,ye)  =  f  g((-xe + yc)/2,  (xe + y^/2)  For  xg  Mod  2  =  ye  Mod  2 

l  undefined  otherwise 


This  mapping  is  illustrated  by  Figure  7.  This  figure  shows  the  correspondence  between  points  in  the 
mapping.  The  dashes  indicate  the  points  which  are  not  defined  in  the  new  filter. 

The  algorithm  for  cascaded  filtering  with  expansion  employs  recursive  application  of  the  V2  expansion 
operation.  Each  expansion  enlarges  the  smallest  distance  between  samples  by  V2  and  alternates  the 
direction  of  that  smallest  distance  between  ±45°  and  0o.90°.  For  this,  we  can  define  a  more  general 
expansion  operator:  F.y^A:{.}.  This  more  general  operator  expands  the  filter  to  the  same  grid  as  an  image 
which  has  been  V2  sampled  k  times. 

Each  application  of  the  V2  expansion  operation  rotates  the  filter  by_45°.  For  a  circularly  symmetric  filter 
this  rotation  has  no  effect  and  we  can  express  an  expansion  of  V2 k  as  k  recursive  applications  of  the 
y/2  expansion. 
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.(1,1)  .(0,1)  (1,1) 

.(-1,0)  .(0.0)  .(1,0) 

.(-1,-1).(0,-1).(1,-1) 

maps  into 

•d.l) 

.(0,1)  -  .(1,0) 

.(-1,1)  -  .(0,0)  -  .(1,-1) 

.(-1,0)  -  .(0,-1) 

(-1.-1) 

Figure  7:  Example  of  mapping  given  by  Ey^{} 

The  general  V2  expansion  operation,  EyjA:{g(x,y)},  may  be  expressed  informally  as  follows.  For  each 
point  (x.y)  at  which  the  filter  gfc  l(x,y)  is  defined,  define  a  new  point  in  g^x.y)  at  (x-y,  x  +  y)  and  copy  the 
value  from  gw(x,y)  into  the  point 

This  mapping  may  be  expressed  more  formally  as  follows:  When  k  is  odd,  the  filter  is  mapped  onto  a  grid 
whose  axes  are  ±45°,  and  whose  smallest  distance  between  samples  is  2k/1.  The  points  on  this  grid  are  those  at 
which 

x  mod  2**+1*/2=  y  mod  2^k+  ^/2  =  0. 
e  e 

For  even  k.  the  expanded  filter  will  be  mapped  onto  a  grid  whose  axes  are  at  0°  and  90°.  The  distance 
between  samples  along  these  axes  will  also  be  2kn.  The  mapping  Ey^-A:{g(x,y)}  may  be  defined  as: 

For  even  k : 

EVT  *te(*.y)}  =  g^Ve) =  |  for  K  mod  2> =  0  and  (ye  mod  2>  =  0 

undefined  otherwise 

For  odd  k : 

E^gCx.,)}  =  ge(x^  =  r  g(^±iil^)fbrx,mod2<‘+1'«  =  y.mod***1*'2 

undefined  Otherwise 


3.2.4  Frequency  Domain  Effects  of  V2  Expansion 


The  V2  expansion  operator  has  a  well  defined  effect  on  the  transfer  function  ofjts  argument.  As  with 
V2  sampling,  a  new  Nyquist  boundary  is  created  which  is  a  45°  rotation  and  a  V2  shrinking  of  the  old 
boundary.  Inside  this  new  Nyquist  boundary  is  a  copy  of  the  old  transfer  function  scaled  smaller  in  size  by  a 
factor  of  V2 .  Outside  this  new  Nyquist  boundary  is  a  reflection  of  die  sealed  transfer  function,  'lliis  is 
illustrated  by  Figure  8  below,  which  shows  the  3  dB  contour  of  a  low-pass  filter  before  and  after  the  expansion 
operation.  Figures  9  and  10  show  plots  of  the  transfer  functions  of  the  Gaussian  low-pass  filter  (R=4,  a =4), 
before  and  after  the  expansion  operation.  Note  the  four  lobes  in  the  comers  of  Figure  10.  These  arc  the 
reflections  of  the  pass  region.  If  these  were  to  show  up  in  the  composite  filter  they  could  cause  a  large 
stop-band  response,  which  would  alter  the  locations  of  peaks  and  ridges  in  the  resulting  band-pass  images. 


b 


Figure  8:  Effect  on  T ransfer  Function  of  E^/j  Expansion 

Operator 


E\/2  {•}  sca*cs  s*ze  of  the  transfer  function  by  V2  so  that  it  approximates  the  larger  Gaussian  filter, 
gfx,y)  within  the  new  smaller  Nyquist  boundary.  That  is 

KeV2 te°(*.y)}}  ~  Hgfay)} 

within 

ir  <  |  u  +  v  |  <  n  The  new  Nyquist  boundary. 

Where  7{}  is  the  transfer  function  [17]. 

For  the  parameter  values  R=4,  a =4  the  pass-band  is  well  within  this  new  Nyquist  boundary,  and  the 
reflection  of  the  pass-band  falls  into  the  stop-band  of  the  previous  filter.  That  is,  outside  of  the  new  Nyquist 
boundary, 


Figure  9:  Transfer  Function  G0(u,v)  for  R  =  4.0,  a  =  4.0  Before  v  2  Expansion 


?{go(x<y)  *  go(x.y)} 

will  be  very  small  (i.e.  less  than  -60  dB  where  the  reflected  nodes  are  present,  for  R=4,  a  =4)  and  thus  the 
product 

?{Ev^{g0(x,y)}}  ^go(x-y)  *g0(x.y)} 

will  also  be  very  small  outside  the  new  Nyquist  boundary.  As  a  result,  the  impulse  response  at  low-pass  level 
2  is  approximated  by 

gjtay)  =  g<,(x.y)  *  go(x,y)  *  Ev^-{g0(x,y)} 

Figure  11  is  a  plot  of  the  transfer  function  of  the  level  2  low-pass  filter.  As  can  be  seen  the  response  in  the 
comers  is  so  small  that  it  does  not  register  in  this  plot.  The  filter  was  constructed  by  convolving  g0(x,y)  with 
itself  (  a  =  4,  R  =  4),  and  then  convolving  an  expanded  version  g„(x,y)  with  this  composite  filter.  Thus  this  is 
the  same  impulse  response  which  would  occur  at  low-pass  level  2  of  a  DOLP  transform  computed  using 
cascaded  convolution  with  expansion.  A  logarithmic  plot  of  the  amplitude  of  G2(u,v)  is  shown  in  Figure  12. 
This  plot  spans  120  db  in  amplitude  with  the  vertical  marked  on  the  left  at  intervals  of  10  db.  The  response  in 
the  comer  regions  arc  attenuated  more  than  100  dB  from  the  peak. 


11:  Filter  G2(u,v)  for  R  =  4.0,  a  =  4.0' 
g2(x.y)  =  8o(x,y)  *  g0(x,y)  *  EN/2{g0(x,y)} 


3.2.5  Complexity  Analysis  of  Cascaded  Convolution  with  Expansion 

Ihc  algorithm  for  cascaded  convolution  with  expansion  is  illustrated  by  the  flow  graph  in  Figure  13.  Its 
computational  complexity  may  be  seen  by  an  analysis  of  the  steps  in  the  algorithm. 

Low-pass  image  0.  Ljx.y)  .is  produced  from  the  original  image,  p(x.y)  ,by  convolution  with  g0(x,y) . 
LJx.y)  =  P(x.y)  *  ge(x.y) 

Band-pass  level  0,  <3>0(x,y),  is  then  produced  by  subtracting  Ljx.y)  from  p(x,y) . 

$0 (x.y)  =  p(x.y)  -  L  Jx,y) 

The  convolution  requires  N  X0  multiplies  and  additions,  and  the  subtraction  requires  an  additional  N 
additions. 

Low-pass  level  1  is  then  formed  by  convolving  low-pass  level  0  with  the  low-pass  filter. 

LJx,y)  =  LJx,y)  *  g0 (x,y) 

Band-pass  level  1  Ls  then  formed  by  subtracting  low-pass  level  1  from  low-pass  level  0. 
e$>l(x.y)  =  L0(x,y)  -  LJx,y) 

As  with  band-pass  level  0,  the  convolution  requires  NX,  multiplies  and  additions  while  the  subtraction 
requires  an  additional  N  additions. 

Low-pass  level  2  is  then  formed  by  convolving  low-pass  level  1  with  an  expanded  version  of  the  low-pass 
filter.  The  expansion  operation  scales  the  filter  larger  by  a  factor  of  s/2  without  increasing  the  number  of 
coefficients. 

Ljx.y)  =  Ljx.y)  *  E J^{gjx.y)} 

Band-pass  level  2  is  then  formed  by  subtracting  low-pass  level  2  from  low-pass  level  1. 

'Ljx.y)  =  LJx.y)  -  LJx,y) 

Since  expansion  does  not  alter  the  number  of  coefficients  this  convolution  also  requires  N  X0  multiplies  and 
additions  and  the  subtraction  requires  an  additional  N  additions. 

Low-pass  level  3  is  then  formed  by  convolving  low-pass  level  2  with  a  twice  expanded  version  of  the 
low-pass  filter.  Two  applications  of  the  expansion  operation  scales  the  filter  larger  by  a  factor  of  2  leaving  the 
original  filter  coefficients  on  a  grid  with  every  other  row  and  column  set  to  zero. 

Ljx.y)  -  LJx,y)  *  E^l{gjx.y)} 

Band-pass  level  3  is  then  formed  by  subtracting  low-pass  level  3  from  low-pass  level  2. 

$  Jx.y)  =  LJx.y)  -  LJx,y) 

Since  expansion  docs  not  alter  the  number  of  coefficients  this  convolution  also  requires  N  X0  multiplies  and 
additions  and  the  subtraction  requires  an  additional  N  additions. 

In  a  similar  manner,  each  band-pass  image  k  is  produced  by  first  creating  low-pass  image  k  by  convolving 
low-pass  image  k-1  with  a  copy  of  the  low-pass  filter  which  has  been  expanded  k-1  times. 


1/x.y)  =  J lk_/x.y)*  E^k-l^g,, (x,y)} 

Low-pass  image  k  is  then  subtracted  from  low-pass  image  k-1  to  produce  band-pass  image  k. 

\(x,y)  =  tk  i(x.y)  -  lk(x.y) 

Since  expansion  docs  not  alter  the  number  of  coefficients  each  convolution  requires  N  X0  multiplies  and 
additions  and  each  subtraction  requires  an  additional  N  additions. 

Since  there  arc  K  =  Log^N)  band-pass  images,  the  total  cost  is 

C  =  X0  N  Log^N)  multiplies  and 
(X0  +  1)N  Logs(N)  additions. 

Since  cascaded  convolution  docs  not  involve  resampling  the  any  of  the  images,  the  memory  costs  for 
computing  a  DOLP  transform  in  this  manner  arc  not  affected.  As  with  equation  (6),  the  memory 
requirements  arc 

M  =  N  l.oggfN)  memory  cells 

3.3  Resampling  and  Cascaded  Convolution  with  Expansion 

The  computational  cost  and  memory  requirements  for  a  DOLP  transform  can  be  reduced  substantially  by 
resampling  each  low-pass  image  before  each  cascaded  convolution.  The  savings  in  computational  complexity 
result  because  there  resampling  reduces  the  number  of  points  at  which  the  convolution  is  evaluated  for  each 
new  level,  while  cascaded  convolution  holds  the  number  of  filter  coefficients  constant  In  this  fast  algorithm 
recursive  expansion  of  the  low-pass  filter  is  not  needed.  In  the  odd  number  levels,  expansion  is  given 
implicitly  by  the  resampling.  In  the  even  numbered  levels,  a'  single  V2  expansion  is  needed  to  place  the  filter 
coefficients  on  the  same  sample  grid  as  the  data. 

3.3.1  The  Algorithm  and  Complexity  Analysis 

The  algorithm  for  resampling  and  cascaded  convolution  with  expansion  is  illustrated  in  the  data  flow  graph 
shown  in  Figure  14.  This  algorithm  runs  as  follows.  Low-pass  and  band-pass  levels  0  and  1  are  computed  as 
described  above  for  cascaded  convolution  with  expansion.  That  is,  low-pass  level  0  is  constructed  by 
convolving  the  picture  with  the  low-pass  filter  gQ(x,y). 

LJx.y)  =  p(x,y)  *  g0(x,y) 

Band-pass  level  0,  <&0(x,y),  is  then  produced  by  subtracting  La(x.y)  from  pfx,y) . 

(x.y)  =  p(x,y)  -  LJx,y) 

Thus  the  band-pass  impulse  response  at  level  0  is 
b0(x,y)  =  8(x,y)  -  gjx,y) 

Low-pass  level  1  is  then  formed  by  convolving  low-pass  level  0  with  the  low-pass  filter. 

Lx(x,y)  =  L„(x.y)  *  gjx,y) 

Band-pass  level  I  is  then  formed  by  subtracting  low-pass  level  1  from  low-pass  level  0. 


Figure  14:  Data  Flow  Graph  foe  Comp< 
Resampling  and  Cascaded  C< 


tt/x -y)  =  LJx,y)  -  L/X.y) 

The  impulse  response  at  band-pass  level  1  is 
b/x,y)  --  So^;  -  (  *gjx.y)) 

Both  band-pass  level  0  and  band-pass  level  1  require  X8  N  multiplies  and  (X0  +  1)  N  additions.  They  each 
produce  N  band-pass  samples. 

For  each  band-pass  level  2  through  K-l,  the  low-pass  image  k-l  is  first  resampled  at  V2  by  the  operation 
Sy^-U.  This  resampling  reduces  the  number  of  sample  points  by  a  factor  of  2  from  the  low-pass  image  at 
k-l.  For  odd  levels,  resampling  leaves  the  data  on  a  cartesian  grid,  and  thus  no  expansion  is  necessary.  The 
low-pass  image  or  level  k  is  thus  formed  by  simply  convolving  the  filter  with  the  low-pass  image  from  level 
k-l. 

LJx.y)  =  1^/x.yJ  *  g0(x.y) 

On  even  levels,  resampling  places  the  data  onto  a  \/2  sample  grid.  To  convolve  an  image  on  a  V2  sample 
grid;  the  low-pass  filter  coefficients  must  be  remapped  to  a  V2  grid  by  the  expansion  operation. 

T =  ^./xyj  *Ey/j{g0('x.yJ} 

In  both  eases  the  band-pass  image  is  then  formed  by  subtracting  the  result  of  the  convolution  from  the 
previous  low-pass  image. 

\(x,y)  =l]l  l(x,y)  -  L^(x,y) 

For  S2  =  V2 ,  each  resampling  reduces  the  number  of  sample  points  by  2,  and  thus  reduces  the  number  of 
multiplies  and  additions  by  a  factor  of  2.  Thus  the  total  number  of  multiplies  and  additions  is  given  by 

C  =  X„  N  ( 1  +  1  4  1/2  +  1/4  4-1/8  +  ...) 

=  3  N  X8  multiplies 
and 

3  N(  X8  +  1 )  additions. 

As  with  the  resampling  algorithm  described  above,  the  total  number  of  memory  cells  required  is 
M  =  3  N 

3.3.2  The  Impulse  Responses  for  Cascaded  Convolution  with  Expansion  and  Resampling 

In  the  cascaded  filtering  algorithms  described  above,  the  band-pass  images  arc  formed  by  subtracting 
adjacent  low-pass  images.  The  band-pass  impulse  responses  arc  thus  equal  to  a  difference  of  low-pass  impulse 
responses  which  arc  produced  by  cascaded  filtering.  Because  a  finite  impulse  response  Gaussian  filter  is  only 
an  approximation  of  the  Gaussian  function,  the  low-pass  impulse  responses  for  levels  1  through  K  are  only 
approximations  of  scaled  copies  of  the  level  0  low-pass  impulse  response. 

The  low-pass  impulse  response  at  level  1  is 
g/x,y)  =  g„(x,y)  *  gjx,y) 

Thus  at  low-pass  level  1,  a  V2  scaling  in  size  of  g8(x.y)  is  approximated  by  the  simple  cascaded  convolution 
ofg8(x,y). 
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Low-pass  level  2  is  formed  by  resampling  low-pass  level  1  at  a  sample  distance  of  V2  and  then  convolving 
with  an  expanded  version  of  the  low-pass  filter  g0(x,y). 

g/x.y)  =  Ey/jlSofay)}  *$y/2{go(x'y)*  zJx.y)} 

The  low-pass  image  from  level  2  is  then  resampled  at  a  distance  of  s/2  for  a  second  time,  which  places  it  on 
*  a  sample  grid  with  a  unit  distance  of  2.  TTiis  low  pass  image  is  then  convolved  with  the  low  pass  filter  gc(x,y). 
The  resampling  provides  a  remapping  of  the  filter  coefficients  and  so  no  expansion  is  needed  at  this  level. 
Thus  the  size  scaling  of  g0  by  a  factor  oils/ 2  is  approximated  by 

g/x.y)  =  go(x-y)  * (F'VT {g0(x.y)}  *s ■yfi{gjx.y)  *  g0(x.y)}} 

In  general,  the  impulse  response  at  low-pass  level  k,  from  k=2  to  K-i  is  given  by  the  following  recursive 
relationships  depending  on  whether  k  is  even  or  odd: 

For  even  k: 

gk(x,y)  =  E^{g0(x.y)}  *Ss/r^{gk_/x,y)} 

And  for  odd  k: 

g/x,y)  =  ga(x,y)  *  S yfi{gk.fxj)} 

3.3.3  The  Size  of  the  Impulse  Responses 

Size  scaling  the  kernel  low-pass  impulse  response  by  resampling  the  continuous  Gaussian  function  at  a 
denser  sample  rate  would  yield  a  sequence  of  radii  Rk  given  by 

Rk  =  R02(k/2) 

The  sequence  of  radii  is  somewhat  different  with  cascaded  filtering.  In  this  ease,  the  expansion  operation 
maps  the  furthest  coefficient,  at  say,  (R,0),  to  a  new  point  at  (R,R).  This  gives  an  increase  in  radius  of  s/2 . 
Convolution  with  the  composite  low  pass  filter  then  adds  this  new  size  to  that  of  the  composite  filter. 

That  is,  at  level  0  the  radius  is  R0.  At  level  1  the  composite  filter  is  the  auto-convolution  of  g0(x,y),  and  its 
radius  is  thus  2R„-1.  The  level  2  composite  filter  is  formed  by  convolving  the  level  1  composite  filter  with  an 
s/2  expanded  version  of  g„.  The  radius  of  the  level  2  composite  filter  is  thus  2R0  +  s/l  R0  -  2.  A  general 
formula  for  the  radius  at  any  level  k  >  0  is 

(k-l) 

Rk  =  R„  -  k  +  R0£  (V2)"'1 

n=Q 


4  An  Example  of  the  DOLP  Transform 

Figure  15  shows  a  resampled  DOLP  transform  of  an  image  of  a  teapot  that  was  produced  using  the  fast 
computation  techniques.  In  this  Figure  the  image  at  the  lower  right  is  the  high  frequency  image,  &9(x,y). 
The  upper  left  comer  shows  the  level  1  band-pass  image,  “BY x.y),  while  the  upper  right  hand  corner  contains 
the  level  2  band-pass  image,  c$>i(x,y).  Underneath  the  level  1  Band-pass  image  are  levels  3  and  4,  then  5  and  6, 


Figure  16:  Levels  5  Though  13  of  the  Resampled  DOLP  Transform  of  a  Teapot  Image 


etc.  Figure  16  shows  an  enlarged  view  of  band-pass  levels  5  through  13.  This  enlargement  illustrates  the 
unique  peaks  in  the  low  frequency  images  that  occur  for  each  grayscale  form. 

These  images  were  formed  using  both  resampling  and  cascaded  convolution  with  expansion.  Each  band¬ 
pass  impulse  response  is  composed  of  a  difference  of  Gaussian  low-pass  filters  with  a  ratio  of  scales  of  S2  = 
V2 .  These  band-pass  images  were  computed  by  forming  low-pass  images  with  the  cascaded  convolution 
with  expansion  technique  and  then  subtracting  to  form  the  Band-pass  images.  'ITic  use  of  V2  resampling  is 
apparent  from  the  reduction  in  size  for  each  image  from  level  3  to  13.  In  the  even  numbered  images,  on  the 
right  of  each  pair,  the  image  is  actually  on  a  V2  sample  grid.  I'o  display  these  V2  images,  each  pixel  was 
printed  twice,  creating  the  interlocking  brick  texture  evident  in  Figure  16. 

5  Summary  and  Conclusions 

This  paper  has  defined  the  Difference  of  Low-Pass  (  DOLP  )  transform.  The  DOLP  transform  is  a 
reversible  transform  that  separates  a  signal  into  a  set  of  band-pass  components.  The  DOLP  transform  serves 
as  the  basis  for  a  representation  for  two-dimensional  shape  that  is  described  in  a  companion  paper  [11].  The 
DOLP  transform  is  shown  to  require  0(N2)  multiplies  and  produce  0(N  Log(N))  samples. 

The  DOLP  transform  is  interesting  because  shapes  (  and  signals )  which  arc  represented  by  encoding  peaks 
and  ridges  (  or  zero-crossings  )  in  the  DOLP  transform  can  be  matched  efficiently  despite  changes  in  size, 
orientation,  or  position,  and  despite  corruption  by  image  noise.  One  of  the  biggest  obstacles  to  use  of  the 
DOLP  transform  for  describing  and  matching  shapes  in  images  was  the  apparent  computational  and  memory 
costs.  In  this  paper  we  have  described  two  independent  techniques  which  may  be  used  to  reduce  the 
computational  complexity  and  storage  costs  of  a  DOLP  transform.  The  technique  of  resampling  is  shown  to 
reduce  the  computational  complexity  of  a  DOLP  transform  to  0(N  Log(N))  multiplies  and  the  storage 
requirements  to  0(N)  samples.  The  technique  of  cascaded  convolution  with  expansion  is  also  shown  to 
reduce  the  computational  cost  of  a  DOLP  transform  to  0(N  Log(N))  multiplies,  but  does  not  affect  the 
storage  requirements.  It  is  then  shown  that  these  two  techniques  may  be  combined  to  produce  a  DOLP 
transform  in  O(N)  multiplies  that  requires  O(N)  samples. 

Cascaded  convolution  has  been  investigated  recently  as  a  technique  for  efficiently  realizing  large  digital 
FIR  filters  [1].  In  particular,  Burt  [5]  has  employed  a  cascaded  convolution  of  a  kernel  which  is  an 
approximation  to  a  Gaussian  to  obtain  larger  "Gaussian-like"  filters.  Such  ajproccss  requires  a  doubling  in 
the  number  of  convolutions  with  the  fixed  size  kernel  for  each  increase  of  V2  in  filter  size.  Our  use  of  the 
expansion  function,  however,  permits  a  composite  Gaussian  filter  of  size  S  V2  to  be  formed  from  a  composite 
Gaussian  of  size  S  by  one  convolution  of  the  kernel  filter.  This  technique  is  general  and  should  be  of  benefit 
whenever  low-pass  kernel  filters  are  cascaded  to  form  larger  impulse  responses. 

The  scale  factor  of  V2  for  filter  size  results  naturally  from  both  fast  techniques.  In  resampling,  it  occurs 
because  it  is  the  smallest  distance  larger  than  one  between  samples  on  a  cartesian  grid.  It  is  the  smallest jrate  at 
which  a  two-dimensional  discrete  sequence  can  be  resampled  without  interpolation.  The  factor  \/2  also 
occurs  with  cascaded  filtering.  It  is  the  increase  in  size  jcalc  provided  by  convolving  a  Gaussian  low-pass  filter 
with  itself.  This  happy  coincidence  indicates  that  V2  is  a  very  convenient  value  for  the  scale  factor  for  a 
DOLP  transform  that  is  to  be  used  to  represent  images  for  matching;  And,  indeed,  this  factor  turns  out  to 
work  quite  well  (10)  for  representation  and  matching  with  the  DOLP  transform. 


The  most  important  result  of  this  work  is  that  it  makes  available  the  representational  power  of  the  DOLP 
transform  without  a  prohibitive  cost  in  computation.  For  a  256  by  256  image,  if  the  separable  form  of  the 
Gaussian  filter  is  used,  the  total  cost  of  computation  for  the  16  band-pass  images  is 

C  =  3  x  18  x  2562  =  3.538  million  multiplies 

compared  to 

C  =  18  x  2564  =  77,309.41133  million  multiplies 

without  the  techniques  of  cascaded  convolution  with  expansion  and  resampling.  Thus,  the  calculation  of  a 
DOLP  transform  in  under  a  second  is  made  possible  by  implementing  these  fast  techniques  on  commercially 
available  high-speed  vector  processing  peripherals. 
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